Experimental realisation of Shor's quantum factoring algorithm using qubit recycling 
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Quantum computational algorithms exploit quantum mechanics to solve problems exponentially 
faster than the best classical algorithm^^^. Shor's quantum algorithrrP for fast number factoring is 
a key example and the prime motivator in the international effort to realise a quantum computei^. 
However, due to^the substantial resource requirement, to date, there have been only four small-scale 
demonstrations^ ^. Here we address this resource demand and demonstrate a scalable version of 
Shor's algorithm in which the n qubit control register is replaced by a single qubit that is recycled 
n times: the total number of qubit s is one third of that required in the standard protoco r°"^^ l 
Encoding the work register in higher-dimensional states, we implement a two-photon compiled 
algorithm to factor N — 21. The algorithmic output is distinguishable from noise, in contrast to 
previous demonstrations. These results point to larger-scale implementations of Shor's algorithm 
by harnessing scalable resource reductions applicable to all physical architectures. 



Shor's factoring algorithm consists of a quantum or- 
der finding algorithm, preceded and succeeded by various 
classical routines. While the classical tasks are known 
to be efficient on a classical computer, order finding is 
understood to be intractable classically. However, it is 
known that this part of the algorithm can be performed 
efficiently on a quantum computer. To determine the 
prime factors of an odd integer one chooses a co- 
prime of A/", X. The order r relates x io N according to x'^ 
mod A" = 1, and can be used to obtain the factors, given 
by the greatest common divisor gcd(x^ ± 1, A'). 

The quantum order finding circuit involves two reg- 
isters: a work register and a control register. In the 
standard protocol, the work register performs modular 
arithmetic with m = \\0g2 N] qubits, enough to encode 
the number A^, and the n qubit control register provides 
the algorithmic output, with n bits of precision. 

Measuring the control register in the computational 
basis will yield a result of k2^/r where k is an integer 
between and r — 1, with the value of k occurring proba- 
bilistically. Dividing the result by 2^ gives the first n bits 
oik/r and r may be found with classical processing, us- 
ing the continued fraction algorithm. For large n, and a 
perfectly functioning circuit, the output probability dis- 
tribution of the control register is a series of well defined 
peaks at values of A:2^/r (Fig. ^p)- (See Appendix for 
details.) 

Here we implement an iterative version of the order 
finding algorithm^^^, in which the control register con- 
tains only a single qubit which is recycled n times, using 
a sequence of measurement and feed- forward operations, 
with each step providing an additional bit of precision 
(Fig. [1^). Reducing the number of qubits in quantum 
simulations and quantum chemistry has been achieved 
with recursive phase estimation^^ while ground state 
projections have been demonstrated by exploiting similar 
techniques in NMRP^^ 

The iterative version of the order finding algorithm, 
displayed in Fig.[l^, is closely related to the semi-classical 
picture of the quantum Fourier transfornP^. Rather than 
performing a Fourier transformation across all control 



register qubits simultaneously, it is sufficient to measure 
the coherence between computational basis states of in- 
dividual qubits in the control register (from least signifi- 
cant bit to greatest significant bit) deciding on the phase 
coherence of the next measurement depending on the pre- 
vious result. In fact, the control register need not contain 
more than one qubit at any one time. This register be- 
gins in a product state and all unitaries controlled from 
the (i)th control qubit are performed before those con- 
trolled from the (i+l)th control qubit. This means that 
all operations on the (i)th qubit can be performed before 
the (i+l)th qubit is initialised and a single control qubit 
can be re-used, or recycled, with the state of the work 
register iteratively updated at each stage. 

To guarantee a good level of precision in /c/r, a full 
scale implementation of Shor's algorithm requires n to 
be log{N'^) < n < log{2N'^). So for full scale imple- 
mentations, qubit recycling reduces the total number of 
qubits required from [3 log2 A^] to [log2 A'] + 1; the only 
penalty is a polynomial increase in computation time, 
while the exponential speedup is retained — i.e. it is scal- 
able. In general, saving in qubits can be more than 2/3 
if more control qubits are required, or less than 2/3 in 
smaller proof of principle demonstrations such as this. 

As the number of control qubits, or iterations, n is 
reduced, in general, the precision is reduced and the k 
peaks in the probability distribution become smeared, as 
shown in Fig. [TJd. However, in the special case where 
the order is a power of two and r = 2^, the peaks cor- 
respond exactly to the logical states of p qubits^^ such 
that the output is equivalent to that of an incoherent 
mixture of p qubits (any additional control qubits re- 
main unentangled throughout the algorithm and simply 
encode the trailing zeros in the uniform distribution). 
Factoring A^ = 15 gives either order 2 or order 4 for each 
of its co-primes^ ^; independent verification of entangle- 
ment is therefore requirecf^. For this reason we focus on 
A" = 21 with the co-prime x = 4 to give ordei^^ r = 3. 

To two bits of precision the expected outcomes 00, 
01, 10 and 11 occur with probabilities |, ^, | and 
|, respectively (Fig lip). Quantum interference in the 



2 




1 00 01 10 11 000 oil 101 0000 0101 1011 



Increasing precision 

FIG. 1: The iterative order finding algorithm for factoring 21. a, Measurement of the control qubit after each controlled 
unitary gives the next most significant bit in the output and the outcome is fed forward to the iterated (semi-classical) Fourier 
transform, which applies either the identity operation / or the appropriate phase gate prior to the Hadamard H. b, As 
the number of iterations increases the precision increases, c, For two bits of precision the controlled unitary operations can be 
constructed with this arrangement of controlled-swap gates. 



Fourier transform is constructive for states contributing 
to the 00 term and boosts its probability of observation 
to three times that of the probability for observing the 
10 term, which experiences destructive quantum interfer- 
ence among its contributory states. Decoherence in the 
Fourier transform would degrade the contrast between 
these terms. For states contributing to the 01 and 11 
terms, the Fourier transform imparts phases to equalise 
the probability of observing these outcomes. Therefore, 
the underlying periodic structure apparent in the two 
qubit probability distribution is susceptible to decoher- 
ence, and therefore distinguishable from noise (See Ap- 
pendix for details). 

We implement a scalable iterative quantum algorithm 
with a compiled version of the quantum order finding 
routine where the circuit is constructed for a particu- 
lar factoring case (here = 21 and x = 4) admitting 
an experimentally tractable implementation, as shown in 
Fig. [2^. There are several st eps t o compilations, common 
to previous demonstrations^^*^^. Firstly, for = 21 and 
X = 4, unitaries and their decompositions can be calcu- 
lated explicitly, as can the full state evolution; secondly, 
redundant elements of these unitaries are omitted. Fi- 
nally, and specific to our own demonstration, since only 
3 of the possible 2^ levels of the conventional 5 qubit 
work register are ever accessed, a single qutrit is used for 
the work register, instead of 5 qubits, to realise a hybrid 
qubit-qutrit systenP^ . (See Appendix for details). 

The controlled unitaries that apply the function 
f{x,a^N) = x^ (mod A^) (where a is a control regis- 
ter computational basis state) to the work register may 
be realised with a sequence of controlled swaps. For our 



two qubit control register, the single swap of U'^ is im- 
plemented with a controlled-NOT (CNOT) gate; is 
realised with two swaps, the first of which is a CNOT 
gate, while it is sufficient for the second swap to be un- 
controlled. (See Appendix for details). 

Our scheme therefore requires two consecutive pho- 
tonic CNOT gates — something that has not previously 
been demonstrated — acting on qubit subspaces of the 
qutrit. In our experimental implementation we use the it- 
erative approach with post-selection in place of measure- 
ment and feed-forward; measurement and feed-forward 
operations have been achieved in the context of cluster 
state quantum computing with photonP^. The quan- 
tum order finding circuit of Fig. [2^ was experimentally 
constructed using the optical circuit shown schemati- 
cally in Fig. (2)3. Realising two consecutive CNOT gates 
on two photons was achieved by using an entangle- 
ment driven CNOT gate (eCNOT)^^, followed by a post- 
selected CNOT gate (pCN0T)PH2l (pCNOT gates can- 
not be used in series without the addition of ancilla pho- 
tons). The circuit was constructed with Jamin-Lebedeff 
polarisation interferometers in a calcite beam displacer 
architecture, chosen to provide interferometric stability. 
Photons were generated with a polarisation entangled 
spontaneous parametric down conversion source^^ (See 
Methods for further experimental details). 

The correct algorithmic output from the quantum or- 
der finding circuit for factoring = 21 is confirmed by 
the data shown in Fig.|3^. The two-qubit control register 
output probabilities for 00, 01, 10 and 11 were measured 
and found to have a fidelity of 99 ± 4%^^ with the ideal 
probabilities |, ^, | and |, respectively. The distribu- 
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FIG. 2: Compiled iterative order finding algorithm, a, Idealised compiled circuit diagram. The first controlled unitary is 
implemented as a CNOT acting on work modes and 2; the final controlled unitary uses a CNOT on work modes and 1 
followed by an uncontrolled swap on work modes and 2. P indicates a projection onto the computational 0, preceded by 
either a bit flip gate X, or the identity operation /; H denotes the Hadamard gate. The iterative Fourier transform includes the 
rotation R = \H) {H\ — i\V) {V\ when the first projection is made on to the computational 1 state (See Appendix for details), 
b, Schematic of the experimental circuit (See Methods for details). 



tion on the first qubit, which determines the second bit 
in the total probability distribution, was measured by 
comparing total control register counts, heralded by the 
W{2) detector, for each setting of the first qubit projec- 



tor (stage (2) of the circuit in Fig. 
uniform with fidelity 99 ± 5% (Fig. 



2 3) and found to be 
SJd). 

Since the non-uniformity of the distribution of Fig. [3^ 
arises from the second iteration of the algorithm when 
the first iteration gives output zero, we now focus on 
this case. Analysis of the circuit and algorithm reveal 
a critical dependence on decoherence, particularly in the 
pCNOT and the Fourier transform phase: phase insta- 
bility drives the output toward a uniform probability dis- 
tribution. We confirmed this analysis experimentally, as 
described below. 

With stage 2 of the circuit set so that the first qubit 
is projected onto the computational zero state. Figs. [Sj^ 
and |3]i show 16 probability distributions for the cor- 
rect phase setting {i.e. 0) to implement the semiclassi- 
cal Fourier transform, and 15 incorrect settings of this 
phas^^. When heralded by W(0, 1), the detector that 
does not distinguish between the and 1 states of the 
work qutrit, the control qubit should be in a maximal 
mixture and therefore insensitive to this phase — i.e. give 
a flat line response. In contrast, since mode 2 of the qutrit 
is not involved in the final CNOT gate, the control qubit 
is not entangled with it. Therefore, when heralded by 
the work qutrit being in the 2 state, the control photon 
should be in a pure state and should exhibit maximum 
sensitivity to the phase in the Fourier transform — i.e. a 



unit visibility fringe. Taking into account the relative 
probability amplitudes in the ideal output state, the ideal 
plot of probability distributions would show two sinu- 
soidal curves of full visibility (from detectors C(0)&iy(2) 
and C{1)&:W{2)) each bisected by a flat line (from detec- 
tors C{0)kW{0, 1) and C{l)kW{0, 1)). For the 8 prob- 
ability distributions between two settings of the Fourier 
transform phase, and tt (indicated by black vertical 
lines) we find an average fidelity between this situation 
and our data of 91.6 ± 0.6%. 

The data in Figs. ^ and [3]i can be used to show 
the sensitivity of the circuit to decoherence: Integrat- 
ing the probability distribution over the range to tt of 
the Fourier transform phase simulates the effect of phase 
instability. The red distribution plotted in Fig. [3^ shows 
the near uniform distribution that results from this pro- 
cedure, confirming susceptibility of this circuit to deco- 
herence. 



Methods 

Photon source: Entangled photon pairs, spectrally de- 
generate at 808 nm, were generated in a Type 1 sponta- 
neous parametric down conversion source, with a 404 nm 
CW laser was focused to a 40 jam waist in a pair of crossed 
BiBaOe (Bismuth Triborate) non linear crystals^. Pho- 
tons were spectrally filtered with high transmission inter- 
ference filters of FWHM 3nm, then collected into polar- 
isation maintaining optical fibres (PMFs). PMFs would 
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FIG. 3: Demonstration of order finding. Ideal probability dis- 
tributions are plotted in solid black lines, a, The two-bit out- 
put probability distribution for the order finding algorithm, 
b, Output from the first iteration of the order finding algo- 
rithm, c & d, Probability distributions for 16 phases of the 
Fourier transform, as described in the text, for control qubit 
detectors and 1 respectively, e. The distribution obtained 
from experimental simulation of decoherence in the Fourier 
transform. Data were corrected for the measured difference 
in total coupling efficiency between the work register detec- 
tors. 



normally decohere the polarisation of photons that are 
not aligned with the slow or fast axis of the fibre, as is 
the case for photons that are entangled. Our fibres were 
cut at the midpoint and spliced together with a 90 de- 
gree twist such that the slow axis in the first length was 
aligned with the fast axis in the 2nd length. While this 
modified fibre imparts an unknown phase shift between 
the two polarisations, their coherence is preserved. The 
unknown unitary is pre-compensated with wave plates 
in the source. State tomography of the photon source, 
which drives the eCNOT gate, revealed a highly entan- 
gled state with fidelity 96.9 ± 0.2% to the corresponding 



Bell state"^. 

Optical Circuit: The optical circuit of Fig. was 
experimentally constructed using an architecture of cal- 
cite beam displacers (BD), which separate the ordinary 
and extraordinary polarisations and can be used to form 
very stable Jamin-Lebedeff interferometers — polarisation 
interferometers with parallel light-paths that provide in- 
terferometric stability. The PBSs and Extended PBSs in 
Fig [2] were directly implemented with a single BD. The 
eCNOT gate^^ requires two non polarising 50% reflectiv- 
ity beam splitters (BS in Fig|2| the unitary operation of 
which was constructed with four BDs and wave plates. 
The action of polarisers Al and A2 was realised using 
beam stops after the first BD in the BSs. The operation 
of the circuit is as follows: 

Stage 1 in Fig. [^: The experimental control and 
work registers are initialised within the eCNOT gate 
by respectively configuring the AO and B{) wave plates 
to output the desired states as if each of their inputs 
were the computational |0): AO is set to implement 
the Hadamard and B{) implements the Identity opera- 
tion. The eCNOT gate is driven by pre-entanglementP^ 
from the polarisation entangled SPDC source in the state 
+ |ly,L) (where H/V denotes hori- 

zontally/vertically polarised light and U/L denotes up- 
per/lower path) which is then converted to path entan- 
glement with polarisation beam splitters (PBS) and po- 
larisation flips (X). After combing the two double-rails 
on non-polarising beam splitters (BS) and post select- 
ing on the cases where photons emerge in the two lower 
paths, the 2x2 transition matrices of the optical elements 
{Al, A2, Bl, B2} combine as Al(g)51+A2(g)52: choosing 
a vertical polariser for Al, a horizontal polariser for A2, 
the Identity operation for 51, and a polarisation flip for 
52, implements the CNOT gate logic on the initialised 
states. In its general form, the eCNOT gate can perform 
any controlled unitary operation (by choosing appropri- 
ate optical elements for A1,A2,51, and 52), and the 
addition of a KLM-like teleportation scheme^ allows the 
gate to work with non separable states. 

The polarisation modes within the control spatial 
mode correspond to the qubit computational states in- 
dicated by C(0, 1); at this point the polarisation modes 
within the work spatial mode correspond to the |0) and 
|2) qutrit states, indicated by VK(0,2). 

Stage 2 in Fig. The control qubit is projected onto 
one of the computational states, dependent upon whether 
/ or X is performed before the upper PBS. The lower 
PBS introduces the third mode for the work register so 
that the |0) and |1) states are polarisation encoded in 
the upper spatial mode of the work register (though at 
this stage the |1) state has zero probability amplitude, 
i.e. vacuum) while the lower spatial mode contains only 
one polarisation and corresponds to the |2) state. 

Stage 3 in Fig. \^'. The pCNOT gate relies on pho- 
tonic quantum interference tuned by the half wave plate 
T which is set to 62.5°. Successful operation is heralded 
when one photon is present in the control modes and one 
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photon is present in the work modes. Here, further bal- 
ancing loss is introduced into the W{2) mode. The out- 
put from W{2) and the usual pCNOT work loss mode 
share the same spatial mode but different polarisations. 
The entangling capability of the pCNOT gate was tested 
with a Bell inequality violation (while in situ) recording 
a CHSH value of 2.67 ±0.01 (violating the classical limit 
of 2 by 55 standard deviations). 

Stage 4 ^W- [H^- The control qubit is assigned a 
phase according to the projector in the first iteration, 



allowing implementation of the semi-classical Fourier 
transform. The control qubit states are individually pro- 
jected and provide the order finding results. At the fi- 
nal stage, the work qutrit plays no role in providing or- 
der finding information (other than to herald the control 
qubit) so individual computational states may be traced 
out in detection. The polarisations of the upper work 
spatial mode are not distinguished, but the remaining 
work mode is; these two cases provide a useful method 
to confirm correct circuit operation. 
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Appendix A: Supplementary Information 

Standard protocol operation: The quantum order finding circuit involves two registers: a work register and 
a control register. In the standard protocol, the work register performs modular arithmetic with m = \\0g2 N] 
qubits, enough to encode the number A^, and the n qubit control register provides the algorithmic output, with 
n bits of precision. The control register, |c), starts in a superposition over all logical states, \cini) = Xla=o^ 1^)- 
The states of the work register, \w)^ are then computed as a function of those in the control register the 
coprime x, and the number to be factored A/", 1^;^) = |/(x,a. A/")) = (mod A")), to produce a highly entangled 
state, \c,w) = Xl^^o 1^) 1^^ (mod N)). The condition = 1 induces a periodicity in the work register so that the 
logical states of the control register can be grouped into r sets as Xlj=o(Xli=o K^ + j)) \^^) (states are not 

normalised here and throughout). Crucially, the periodicity in control register states is exploited through a Fourier 
transformation, which transforms a function with period r into a new function with period 2'^/r. 

Compiling the algorithm: The circuit schematic shown in Fig. 2a is designed to perform factoring for the specific 
case of A' = 21 and x = 4. Controlled unitary operations perform f{x^a^N) on the work register for each value of a 
in the control register; for A^ = 21 and x = 4, the unitaries, their decompositions, and the full state evolution, can 
be calculated explicitly by hand and redundant elements of these unitaries are omitted (see next section for further 
details). In the standard qubit encoding of Fig. la, the work register requires 5 qubits to represent A^ = 21 in binary. 
However, only three (of the possible 2^) work register states are ever accessed since the function /(x,a, A^) = x^ 
(mod A^) is a periodic repetition of {4^ = 1,4^ = 4,4^ = 16}. We implement the work register using a single qutrit, 
taking advantage of the further degrees of freedom available in photon path and polarisation encoding. The qutrit 
levels represent the active sates {1,4, 16} with the qutrit labels taking the log4, of these values. This replaces the 5 
qubits shown in Fig. la. 

Unitary decomposition and algorithm operation: The function /(x, a, N) = x^ (mod N) is carried out on the 
work register through a sequence of controlled unitary transformations {U^^ ^}, with each unitary controlled from 
the (j)th control qubit (where j = 1 is the final qubit, giving the most significant bit). With the the work register 
initialised in the \wini) = |1) state, should perform the mapping (mod A"). Therefore, f{x^a^N) 

is realised through a sequence of controlled permutations, and each controlled permutation may be realised with a 
sequence of controlled swaps. For A^ = 21 and x = 4 with two control bits, U'^ : {1 ^ 16,4 ^ 1,16 ^ 4} is 
implemented first, followed by : {1 ^ 4, 4 ^ 16, 16 1}. It is shown in Fig. 2a that, since the standard work 
register starts in the state |1), only the mapping |1) |16) from is required; with our relabelling for the work 
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qutrit, the work register initialised as |0) should be mapped to |2), which is realised with a CNOT gate acting on a 
qubit subspace of the qutrit: |0) O |2). may be performed with two swaps: 0^1 followed by ^ 2. The first of 
these swaps is implemented as a CNOT, while an uncontrolled swap, equivalent to a relabelling of two modes of the 
qutrit, is sufficient for the second (see below for further details). 

The (still unnormalised) state of the qubit-qutrit system shown in Fig. 2a before the first projection is \ci,^Wt)i = 
|0) (|0) + |2)) + |1) (|0) — |2)), so that selecting either the identity operation / or the bit flip X before the polariser P 
post selects on to the first or second term respectively. Ideally, there is an equal probability of observing |0) or |1) in 
the control register, which corresponds to the least significant digit in the final output. With |0) post selected in the 



control register, the state after the second CNOT is 



.(0) 



|00) + 111) + 102) + 112); the state after the subsequent 



action of a controlled swap cS on work register modes and 2 is found to be equivalent to that after an uncontrolled 
swap uS on the same modes: cSo^2 '^t / = I ^ uSo^2 '^t / = |00) + |10) + |11) + |02). When the first qubit 



.(0) 

is projected onto |0), the second qubit component of the Fourier transform is applied by the Identity operation / 



followed by the Hadamard so that the final state is 



.(0) 



fin 



|00) + |00) + |01) + |02) + |10) - |10) - 111) + |12). 



The crucial role of quantum interference in the Fourier transform can now be seen: the probability of observing 
the term in the control register is boosted by constructive interference and is three times that of the probability of 
observing the 1 term, which suffers destructive interference; these terms are the digits 00 and 10 in the final probability 
distribution. Therefore, contrast between these terms degrades with decoherence in the Fourier transform. 

When the ffist control qubit has been projected onto |1), the state after the second CNOT is c'^\w^ = |00) + 

1 11) — 1 02) — 1 12). The subsequent action of a controlled swap is equivalent to that of an uncontrolled swap up to a 



phase flip on terms with a 1 in the control register: / uSq^2 



.(1) 



|00) - 1 11) - |02) + |10). The phase ffip is 



undone at the second qubit component the Fourier transform by applying —i rather than i so that the final state is 
,wt) = |00) - i |00) + i |01) - |02) + |10) + i |10) - i |11) - |12). Here, the Fourier transform imparts phases 



.(1) 



fin 



such that there is an equal probability of observing or 1 in the control register. These terms are the digits 01 and 
11 in the final probability distribution. ^ 



